/*! 
 ***************************************************************************
 * \file   hoffman.cpp
 * Another example of Laplace equation in 2D.
 ***************************************************************************
 *  \author    Luis M. de la Cruz Salas [ Fri Mar 21 12:26:33 GMT 2010 ]
 ***************************************************************************
 * This code solves the Laplace equation in two dimensions: 
 * \f$ \nabla^2 T = 0\f$  in the rectangle \f$ [0,1] \times [0,2] \f$. 
 * The boundary conditions are as shown in the next figure:
 * \image html LaplaceSquare.png "Unit square in 2D" width=10cm
 * \image latex LaplaceSquare.eps "Unit square in 2D" width=10cm
 * \par
 * The exact solution of this problem is: 
 * \f[ 
 * T(x,y) = 10\sin(\pi x) \frac{e^{\pi x}-e^{-\pi y}}{e^{\pi}-e^{-\pi}}
 * \f]
 * \par
 *  To compile and run this example type the next commands:
 * \par
   \verbatim
   % make
   % ./hoffman \endverbatim
 * 
 *  The results of this example are stored in \c DataHoffman and are plotted
 *  using OpenDX with the script \c visdat.net, look at the next figures.
 * \image html  HoffmanExa_01.png "Exact solution" width=10cm 
 * \image html  HoffmanNum_01.png "Numerical solution" width=10cm 
 * \image html  HoffmanErr_01.png "Global error" width=10cm 
 */

#include <iostream>
#include <cmath>
#include "Meshes/StructuredMesh.hpp"
#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"
#include "Schemes/LaplaceScheme.hpp"
#include "Solvers/TDMA.hpp"
#include "Utils/num_utils.hpp"
#include "Utils/inout.hpp" 
#include "Utils/GNUplot.hpp"
using namespace Tuna;

// Scalar field in two dimension
typedef TunaArray<double, 2>::huge ScalarField2D;

void exactSolution(ScalarField2D &T, double hx, double hy, double A,
		   double dx, double dy)
{
  double x, y, c1, s1, s2, s3, sum, exact;
  int Nmax = 20;
  int NI = T.ubound(firstDim);
  int NJ = T.ubound(secondDim);
  
  for(int i = 1; i < NI; ++i)
    for(int j = 1; j < NJ; ++j) {
      sum = 0.0;
      for(int n = 1; n < Nmax; n += 2) {
	x = i * dx;
	y = j * dy;
	c1 = 2.0 / (n * PI);
		s1 = sinh( n * PI * (hy - y) / hx );
		s2 = sinh( n * PI * hy / hx);
		s3 = sin( n * PI * x / hx);
		sum += c1 * s1 * s3 / s2;
      }
      T(i,j) = 2.0 * A * sum ;
    }
}

int main() 
{
  double length_x, length_y, Amp;
  int num_nodes_x, num_nodes_y;
  double tolerance;    
  int tdma_iter;
  timer time;
  
// Reading data from the input file
  std::ifstream input_cin ("input_hoffman"); 
  input_cin >> length_x;
  input_cin >> length_y;
  input_cin >> num_nodes_x; 
  input_cin >> num_nodes_y; 
  input_cin >> Amp; 
  input_cin >> tdma_iter;
  input_cin >> tolerance;
  input_cin.close();

  std::cout.precision(10);
  std::cout.setf(ios_base::fixed,ios_base::floatfield);

  StructuredMesh<Uniform<double, 2> > mesh(length_x, num_nodes_x,
				           length_y, num_nodes_y);
  double dx = mesh.getDelta(X);
  double dy = mesh.getDelta(Y);
  mesh.print();

  ScalarField2D T (mesh.getExtentVolumes());
  ScalarField2D Ta(mesh.getExtentVolumes());
  ScalarField2D Tn(mesh.getExtentNodes());
  
  Range all = Range::all();
  Tn(all, Tn.lbound(secondDim)) = Amp;
  InOut::writeToFile(Tn, 0, "./DataHoffman/temp.", dx, dy);
  NumUtils::interpolateToVolumes(T, Tn);
  
  exactSolution(Tn, length_x, length_y, Amp, dx, dy);
  InOut::writeToFile(Tn, 0, "./DataHoffman/exac.", dx, dy);
  NumUtils::interpolateToVolumes(Ta, Tn);
 
#ifdef WITH_GNUPLOT
    GNUplot plotter;
    plotter("set grid");
    plotter("set contour");
    plotter("set view 60,45,1,1");
    plotter("set ticslevel 1.5");
    plotter("set cntrparam level discrete 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95");
    plotter("splot [0:1][0:2][0:100] \"./DataHoffman/exac.00000\" title \"Exact solution\" w l");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    int pause; std::cin >> pause;
    plotter("splot [0:1][0:2][0:100] \"./DataHoffman/temp.00000\" title \"Initial conditions\"w l");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    std::cin >> pause;
#endif

    SparseMatrix< Diagonal<double, 2> > A(num_nodes_x, num_nodes_y); 

    ScalarField2D             b(num_nodes_x, num_nodes_y);
  
    ScalarEquation< LaplaceScheme<double, 2> > laplace(T, A, b,
						       mesh.getDeltas());
    laplace.setDirichlet(TOP_WALL);
    laplace.setDirichlet(BOTTOM_WALL);
    laplace.setDirichlet(LEFT_WALL);
    laplace.setDirichlet(RIGHT_WALL);
    laplace.setGamma(1.0);
    laplace.print();

    time.tic();        
    laplace.calcCoefficients();
    int iter = Solver::TDMA2DY(laplace, tolerance, tdma_iter);
    std::cout << "\n Elapsed time for solving Ax = b : " 
	      << time.toc() << " seconds \n\n";   
    
    std::cout << "\n Iterations = " << iter
	      << "\t Residual = " << laplace.getResidual();

    laplace.update();  
    NumUtils::interpolateToNodes(Tn, T);
    InOut::writeToFile(Tn, 1, "./DataHoffman/temp.", dx, dy);

    cout << "\n\n Error L1  = " << NumUtils::calcErrorL1( T, Ta )
	 << "\t Error RMS = " << NumUtils::calcErrorL2( T, Ta );

    double ReL1, ReL2, ReLmax;
    double AeL1, AeL2, AeLmax;
    double NL1, NL2, NLmax;
    NL1 = NumUtils::normL1(Ta);
    NL2 = NumUtils::normL2(Ta);
    NLmax = NumUtils::normLmax(Ta);

    AeL1 = NumUtils::calcAbsoluteError_L1( T, Ta);
    AeL2 = NumUtils::calcAbsoluteError_L2( T, Ta);
    AeLmax = NumUtils::calcAbsoluteError_Lmax( T, Ta);
    ReL1 = AeL1 / NL1;
    ReL2 = AeL2 / NL2;
    ReLmax = AeLmax / NLmax;
    
    cout << "\n Error L1  Absolute = " << AeL1
	 << "\t Relative = " << ReL1
	 << "\t Norm = " << NL1;

    cout << "\n Error L2  Absolute = " << AeL2
	 << "\t Relative = " << ReL2
	 << "\t Norm = " << NL2;

    cout << "\n Error Lmax  Absolute = " << AeLmax
	 << "\t Relative = " << ReLmax
	 << "\t Norm = " << NLmax;

    for(int i = 0; i <= T.ubound(firstDim); ++i)
      for(int j = 0; j <= T.ubound(secondDim); ++j) {
	T(i,j) = abs(T(i,j) - Ta(i,j));
      }

    
 
    NumUtils::interpolateToNodes(Tn, T);           
    InOut::writeToFile(Tn, 0, "./DataHoffman/erro.", dx, dy);
    
#ifdef WITH_GNUPLOT
    plotter("splot [0:1][0:2][0:100] \"./DataHoffman/temp.00001\" title \"Numerical Solution\" w l");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    std::cin >> pause;
    plotter("splot [0:1][0:2][0:100] \"./DataHoffman/erro.00000\" title \"Error\" w l");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    std::cin >> pause;
#endif
    
    return 0;
}




